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7-H , The quasilinearization method (QLM) of solving nonlinear differential equations is applied to the quantum mechanics by 
casting the Schrodinger equation in the nonlinear Riccati form. The method, whose mathematical basis in physics was discussed 
recently by one of the present authors (VBM), approaches the solution of a nonlinear differential equation by approximating 
the nonlinear terms by a sequence of the linear ones, and is not based on the existence of some kind of a small parameter. It 
q . is shown that the quasilinearization method gives excellent results when applied to computation of ground and excited bound 
state energies and wave functions for a variety of the potentials in quantum mechanics most of which are not treatable with the 
help of the perturbation theory or the 1/N expansion scheme. The convergence of the QLM expansion of both energies and 
wave functions for all states is very fast and already the first few iterations yield extremely precise results. The precison of the 
wave function is typically only one digit inferior to that of the energy. In addition it is verified that the QLM approximations, 
unlike the asymptotic series in the perturbation theory and the 1/N expansions are not divergent at higher orders. 

PACS numbers: 02.30.Mv, 04.25.Nx, 11.15.Tk 



^ : 

Ph. 

Q-r I. INTRODUCTION 

£: 

Realistic physical calculations usually are impossible without different approximation techniques. Correspondingly 
^ . expansions in small parameters, statistical, variational and majority of numerical methods belong to the arsenal of 
^ ' the modern physics. 

Since many equations of physics are either nonlinear or could be cast in the nonlinear form, the possibility of 
^ ■ adding to this arsenal an additional very powerful approximation technique applicable to nonlinear problems was 
pointed out in a series of recent papers [[[J ||, It is called the quasilinearization method (QLM) and its iterations 
are constructed to yield rapid quadratic convergence and often monotonicity. The quasilinearization method was 
i— i ' developed many years ago in theory of linear programming by Bellman and Kalaba |5| as a generalization of the 
Newton-Raphson method || |?j to solve the systems of nonlinear ordinary and partial differential equations. Its 
modern developments and examples of applications to different fields of science and engineering are given in a recent 
[ monograph 

■ In the original works of Bellman and Kalaba j|, ||, however, the convergence of the method has been proven only 
under rather restrictive conditions of small intervals and bounded, nonsingular forces Q| which generally are not 
fulfilled in physical applications. This could explain an extremely sparse use of the technique in physics, where only 
a few examples of the references to it could be found ||, [l0| 0, [l^, |l3). Recently, however, it was shown [Q by one 
of the present authors (VBM) that a different proof of the convergence can be provided which allows to extend the 
applicability of the method to realistic forces defined on infinite intervals with possible singularities at certain points. 
^ , This proof was generalized and elaborated in the subsequent works [|[ || . 

. ^ ■ In the first paper of the series jlj , the analytic results of the quasilinearization approach were applied to the nonlinear 
Calogero equation |)| for the scattering length in the variable phase approach to quantum mechanics, and the results 

J^"*' were compared with those of the perturbation theory and with the exact solutions. It was shown that the n-th QLM 
approximation sums exactly 2™ — 1 terms of the perturbation theory while a similar number of terms are summed 
• • . approximately. The number of the exactly reproduced perturbation terms thus doubles with each subsequent QLM 

. £h ' approximation, which, of course, is a direct consequence of a quadratic convergence. 

The numerical calculation of higher QLM approximations to solutions of the Calogero equation with different 
singular and nonsingular, attractive and repulsive potentials performed in the next work Q has shown that already 

?3 | the first few iterations provide accurate and numerically stable answers for any values of the coupling constant and 
that the number of iterations necessary to reach a given precision increases only slowly with the coupling strength. It 
was verified that the method provides accurate and stable answers even for super singular potentials for which each 
term of the perturbation theory diverges and the perturbation expansion consequently does not exist. 

In the third paper of the series || the quasilinearization method was applied to other well known typical nonlinear 
ordinary differential equations in physics, such as the Blasius, Duffing, Lane-Emden and Thomas-Fermi equations 
which have been and still are extensively studied in the literature. These equations, unlike the nonlinear Calogero 
equation for the scattering length Q considered in references jj], , contain not only quadratic nonlinear terms but 
various other forms of nonlinearity and not only the first, but also higher derivatives. It was shown that again just a 
small number of the QLM iterations yield fast convergent and uniformly excellent and stable numerical results. 
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The goal of the present work is to apply the quasilinearization method to quantum mechanics by casting the 
Schrodinger equation in the nonlinear Riccati form and calculating the QLM approximations to bound state energies 
and wave functions for a variety of potentials, most of which are not treatable with the help of the perturbation 
theory or the 1/N expansion scheme. We show that the convergence of the QLM expansion for both energies and 
wave functions is very fast and that already the first few iterations yield extremely precise results. In addition it is 
verified that the higher QLM approximations, unlike those in 1/N expansion method, are not divergent at any order. 

The paper is arranged as follows: in the second chapter we present the main features of the quasilinearization 
approach to the solution of the Schrodinger equation, while in the third chapter we consider the application of the 
method to computations for the Coulomb, Hulthen, Poschl- Teller, logarithmic, double-well, anharmonic oscillator, 
linear and different power potentials such as r 3 / 2 and r 5 . The results and their comparison with other calculations, 
convergence patterns, numerical stability, advantages of the method and its possible future applications are discussed 
in the final, forth chapter. 



II. QUASILINEARIZATION APPROACH TO THE SOLUTION OF THE SCHRODINGER EQUATION 

The quasilinearization method (QLM) solves a nonlinear n-th order ordinary or partial differential equation in N 
dimensions as a limit of a sequence of linear differential equations. The idea and advantage of the method is based 
on the fact that linear equations can often be solved analytically or numerically using superposition principle while 
there are no useful techniques for obtaining the general solution of a nonlinear equation in terms of a finite set of 
particular solutions. 

The main feartures and equations of the method, appropriate for physics applications, are summed in Refs. 0, ^, 
In this paper we will follow these references since the derivation there is not based, unlike the derivation in Refs. 

|^] , on the smallness of the interval and on the boundedness of the nonlinear term and its functional derivatives, 
the conditions which usually are not fulfilled in physics. 

We would like to use the method in quantum mechanical calculations with the central potential V(r). In order to 
do this we have to rewrite the corresponding radial Schrodinger equation 



2m 



X"(r) 



V(r) 



l{l + i)h 2 

2mr 2 



X(r) = Ex(r) 



(1) 



in nonlinear form. Here x( r ) — vR(r) and R is the radial part of the wave function. Setting h = 1, k 2 = 2m\E\, and 
U(r) = 2mV(r) + 1(1 + l)/r 2 , we obtain the bound state and scattering Schrodinger equations 



d 2 x(r) 
dr 2 



(k 2 + U(r)) X (r) = 0, E <0 



(2) 



and 



d 2 x(r) 
dr 2 



+ (k 2 - U{r)) X (r) =0, E > 



with the boundary conditions at the origin 



(3) 



X(r) ~ r' 

r— >0 



1 + 1 



(4) 



and at the infinity for the potentials falling off at large r 



X (r) ~ e~ Kr , E<0, (5) 

r — >oo 

nl 

X(r) ~ sm(nr - — +6i), E > 0. (6) 

r — >oo Z 

For potentials behaving at large r as A 2 In r/R or X 2 r p with positive R,p and A the boundary conditions at infinity 
should be changed respectively to 
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X(r) ~ e -*f r V^JK dr ( 7 ) 



r — > oc 



I \ _^2A_ p/2+1 

x (r) ~ e . (8) 

r — too 

The boundary condition (|8|) with r changed to |r| holds at both boundaries r = ±oo in the one dimensional 
problem for the double-well potential (r 2 — 16) 2 /128 considered among others in Ref. (m) where we look for both 
ground (symmetric) and first excited (antisymmetric) solutions. 

It is easy to check that the inverse dimensionless logarithmic derivative <f> of the wave function, 

^0) = K ^Aj> ( 9 ) 

satisfies the nonlinear Riccati equations 

= 1 - (1 + W{x))cj) 2 {x), E < (10) 



dx 
and 



#<X< TT', .,,2 



dx 



= 1 + (1 - W(x))r(x), E>0, (11) 



where x = nr and W(x) — U(x/k)/k 2 are the dimensionless variable and potential, respectively. 

To avoid poles of <j>(x) at the bound state energies it is convenient to define (see jl], |) and the references therein) 
a new function u(x) with the help of the equation 

4>{x) = — tanu(x). (12) 

The corresponding equations for u{x) have the forms 

du(x) 
dx 

and 



-1 + (2 + W{x)) sin 2 u{x), E < (13) 



du ( x > . , , it-, ... • - 



dx 



= -1 + W{x) sin 2 u(x), E > 0. (14) 



Similar types of equations were derived earlier by Drukarev |L5| , Bergmann jl6[ , Olsson [JETj , Kynch JT^] , Franchetti 
||, Spruch @, Dashen ||, Calogero § and Babikov ||. 
The boundary conditions for the function u(x), in view of Eqs. (Hfl), respectively, reduce to 



and 



u(x) -=-, (15) 

u(x) ~ ^ -rnr, E < (16) 



u(x)+a; - ^--S l ,E>0. (17) 



The boundary conditions (fa) and rtq) which for the inverse logarithmic derivative 4>(x) have the forms 
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«*>,~.-v£ir°' (18) 

and 

x -p/2 

>0 , (19) 

x — >OG A 

respectively, for the function u(x) therefore read 

u(x) ~ -nir. (20) 

x — >oo 

Here and in Eq. (|l^) n obviously denotes the number of the excited state with n — 1 corresponding to the ground 
state, n — 2 to the first excited state etc. The minus sign in front of n follows from the fact that in the regions of r 
where V(r) < E in view of Eqs. ( p"3| ) and (|^) the derivative du/dx is negative and u(x) is decreasing. Since its value 
at the origin is zero u{x) stays negative which determines the sign in front of n. 

Returning to the variable r and defining a new function a(r) which has the dimension of length with the help of 
the relation (j>(x) — n{r + a(r)) and substituting it into Eqs. ( |i~0| ) and fill ) we obtain the equations 

da ( r ) _ { , 2 i TTf \\i , I ^2 



dr 

and 



= -(k 2 + U(r))(r + a(r)Y, E<0 (21) 



which are very similar to the Calogero equation 

da(r) 
dr 



(k 2 - U{r)){r + a(r)Y, E>0 (22) 



-2mV(r) {r + a(r)f, (23) 



where a{r) has the meaning of the variable s-wave scattering length ||. These equations are obviously a generalization 
of the Calogero Eq. ( |23| ) for arbitrary values of I and k and reduce to it when I and k are equal to zero. 

The QLM prescription jl], ||, ^, ^, || determines the k + 1-th iterative approximation uj < . + i(x) to the solution of the 
first order nonlinear equation in one variable 

^l = f(u(x),x),u(0) = (24) 
dx 

as a solution of the linear equation 

u k+i( x ) = /(■"*, a;) + («fc + i(a;) - u k (x))f u (u k , x), 

«k+i(0) = 0, (25) 

where the functional f u {u,x) — df(u,x)/du is a functional derivative of the functional f(u(x),x). 
The analytical solution of this equation is 



u k+ i(x) = J ds(f(u k (s),s) - f u (uk(s),s)u k (s)) 

xexp / dtf u (u k (t),t). (26) 



The sequence u k (x), k = 0, 1, 2, ... of QLM iterations satisfying Eqs. (|2!j) and (26), converges uniformly and quadrat- 
ically to the solution u(x) of Eq. ( p4| ) if the initial guess for the zeroth iteration is sufficiently good. In addition, 
for strictly convex (concave) functionals f(u(x),x) the difference u k +i(x) —u k {x) is strictly positive (negative) which 
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establishes the monotonicity of the convergence from below (above), respectively. The exact conditions of the con- 
vergence and the monotonicity for the realistic physical conditions of forces defined on infinite intervals with possible 
singularities at certain points are formulated in Ref . (lj . One can also prove that in the quasilinear approximation 
the energy in the Schrodinger equation satisfies the Rayleigh-Ritz variational principle which ensures the quadratic 
convergence in the QLM energy computations. 

We will limit ourselves here to the bound state calculations with Eqs. (|lj) for the negative energy bound states 
and ( |l4| ) for positive energy bound states which are somewhat more complicated than scattering calculations since in 
the former case the boundary condition at infinity determines a discrete spectrum. 

For the negative energies, Eq. (|l3|), the functionals /(u(x), x), F(u(x), x) = / u (u(ir), x) and G(u(x), x)=f(u(x), x) — 
u{x)f u {u{x),x) are given by 

f(u(x),x) = -l + {2 + W(x))sin 2 u(x), (27) 
F{u{x),x) = (2 + W(x))sm2u(x), (28) 

and 

G(u(x),x)= -1 + (2 + W(x)) sinu(x) 

x [sinu(x) — 2u(x) cosu(x)], (29) 

so that Eqs. ( p5|) and ( |26| ) respectively have the forms 

u 'k+i( x ) - u k+ i(x)F(u k {x) 1 x) = G(u k (x),x), (30) 

and 

u k+ i(x) = J dsG(u k (s),s)exp J dtF(u k (t),t). (31) 
For the positive energies the same equations (|27|)-(p9|) hold with (2 + W(x)) replaced everywere by W(x). 



III. QLM BOUND STATE CALCULATIONS AND THEIR COMPARISON WITH THE l/N 
EXPANSION METHOD AND EXACT SOLUTIONS 

In the previous chapter we have cast the Schrodinger equation in the nonlinear Riccati form and wrote the lin- 
ear equations and the boundary conditions appropriate for the bound state calculations with the quasilinearization 
method. 

In this chapter we consider examples of different singular and nonsingular attractive interactions which, in view of 
their large coupling constants, are not treatable with the help of the perturbation theory and for most of which the 
1 /N expansion series are asymptotically divergent as has been shown in Ref. |]l4| . 

Namely, we apply the quasilinearization method to computations with the Coulomb, Hulthen, Poschl- Teller, log- 
arithmic, anharmonic oscillator, linear and different other power potentials such as r 3 / 2 and r 5 as well to the one- 
dimensinal double-well potential (r 2 — 16) 2 /128, and we compare the wave functions and the bound state energies 
obtained by the quasilinearization method (QLM) with their exact values and with results obtained in the l/N ex- 
pansion theory. To show that the method works equally well also for excited states we calculate in the Coulomb, 
linear and double well potentials the first few excited states as well. 

The calculations were done using the differential formulation, Eq. (^5|), of the QLM iteration, for the simple reason 
that the adaptive numerical integration ]23|] together with interpolation proved faster than the integral formulation 
( ^6| ) , mainly due to the processor time taken by the evaluation of the exponential in Eq. ( p6| ) . For each QLM iteration 
number k, k = 0, 1, 2, 3, . . . k m , numerical integration was performed from x = to the matching point x — x m and 
from the upper bound x — x\j to x ~ x m . 

Let us denote the set of iteration-integration parameters by V = {k m , xjj, iVj, . . .}, where k m is the maximum QLM 
iteration index, x\j is the upper bound of the interval, and Ni is the number of interpolation points in each of the two 
subintervals (0,x m ) and (a; m ,a;u). 

The computation was done in two steps. In the first step, x m , the starting values of parameters, Vq, and a k value, 
Ko, near the expected eigenvalue were prescribed. On the last QLM iteration (fc = k m ) the absolute difference between 
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the left-hand side (LHS) and right-hand side (RHS) solutions, Dp (Ko;x m ) = \u k °' (x m ,Ko) — u k ' (x m ,Ko)\, 
was calculated. The whole process was then performed with a new set of parameters V\, where k m , x\j, etc. were 
increased. This was repeated until some number M of steps, when D-p M (kq; x m ) was stabilized to a required accuracy. 

In the second step, the parameter set Vm thus optimized was used to find the zero of Dp> M (n; x m ) — as a function 
of k: the QLM iteration was first performed for two k values lying on opposite sides of the expected eigenvalue, and 
the QLM iteration (k = 1, 2, 3, . . .) was then repeated for each new n value until D-p M {n; x m ) = 0. 

In this process the value of x m was kept constant, which had the consequence that the RHS interval (x m ,x\j) was 
increasing. Both solutions tend to become unstable near x = x m on their respective sides, if the respective interval 
is too large. It turned out however that it was possible to leave x m unchanged, except that x m typically had to 
increase with the number of the excited state. On the other hand, as is evident from the Figures, the starting values 
of parameters (Pq), in particular k m , had to be large enough to overcome the divergent behavior of the solutions near 
x = x m already for the QLM iteration using Vq. It also turned out that the RHS solution quickly assumes the correct 
value, thus allowing reasonably small xtj, and actually making the process rather independent of the exact value used 
for the boundary condition at infinity. 

The precision was controlled in the following way. The differential equation solver |23| was required to return u(x) 
with the precision of the order of 10~ P ° DE . The required precision of Dp during the optimization of V was 10~ Pv . 
Pode w as taken to be larger than Pp by 1 to 3 to test stability. 

The results of the calculations are summarized in Table | and in Figs. |l}|l5|. The calculations are done for the 
s-states since the calculations for p, d states etc. have the same degree of difficulty and could be performed in a 
similar fashion. In the caption of the table V(r) is the potential and n denotes the number of the excited state; m is 
the mass of the particle and is given different values for different potentials in order to enable comparison of the QLM 
bound state energies with those obtained by the 1/N expansion method in Ref. pi] where m takes on values m = 1 
or 1/2 depending on the interaction. In the graphs of the convergence of Uk(x) with iteration index k we present 
for clarity only those iterations which are distinguishable from the final solution (k = k m ) on the graphs; the actual 
number of iterations is higher in order to achieve greater wave function precision. Figures which display the absolute 
differences between successive iterations, |itfc(r) — Uk—i(x)\, or the differences between the successive iterations and 
the exact solution, |iifc(r) — u exac t(a:)|, show the results for the respective optimized parameters sets, D-p M , and for 
the last k value, i.e., at the E of the eigenvalue. 

The required precision of u(x), or the wave function, was Pode = 9 in all cases except in the logaritmic potential 
case where Pode = 6. The number of digits in the values of E in Table || is the number of stable digits when Pp was 
increased up to Pode + 3, except in the cases of Coulomb, Poschl- Teller and Hiilthen potentials, where we display an 
additional (the first incorrect) digit. 



From Table | and Figs. p|]l5| one can conclude that QLM is extremely precise. Energies and the wave function 
for both ground end excited states typically converge to the order of 10 significant digits after about k m = 10 — 20 
iterations though the precision of E is about one digit more than the precision of the wave function. We used 
the numbers of QLM iterations k m such that the precision of the iteration itself, shown by the Figures displaying 
\uk(r) — Uk-i(x)\, was up to 10~ 15 . 



Our calculations confirm numerically the conclusion following from the proof in Ref. |3| that once the quasilincar 
iteration sequence starts to converge, it will continue to do so, unlike the perturbation expansions in powers of the 
coupling constant or in powers of 1/N, which are often given by the asymptotic series and therefore converge only up 
to a certain order and diverge thereafter. In particular, the 1/N expansions of the binding energy of different ground 
and excited states given in Table |, are strongly divergent for logarithmic, double-well, anharmonic oscillator, linear, 
r 3 / 2 and r 5 potentials at orders of about 20 or higher or even before this as it was shown recently by Bjerrum-Bohr 



Based on our results of the QLM computations of the wave functions and bound state energies for many different 
potentials, one can deduce the following important features of the quasilinearization method in the quantum mechanics: 

(i) The quasilinearization method solves the Schrodinger equation by rewriting it in the nonlinear Riccati form and 
by approximating the nonlinear terms by a sequence of the linear ones. It is not based, unlike perturbation or 
1/N expansion theories, on the existence of some kind of small parameter. 

(ii) The quasilinearization method works equally well for both ground and excited states. It is extremely precise: 
binding energies and the wave functions converge to the order of 10 significant digits after about 10—20 iterations. 
Typically, the numerically obtained precison of the wave function is only one digit inferior to that of the energy. 




IV. CONCLUSION 
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(iii) Fast convergence of the QLM iterations to the exact solution confirms numerically the uniform and quadratic 
law of convergence proved in Refs. jl| @, || for realistic physical interactions defined on infinite intervals with 
possible singularities at certain points of the intervals. 

(iv) For convergence it is enough that an initial guess for the zeroth iteration is sufficiently good. In all the examples 
considered in the paper the simplest initial guess of setting the logarithmic derivative of the wave function 4>(x) 
equal to zero or to x at the origin was enough to produce a rapid convergence. 

(v) By using the high numbers of QLM iterations k m such that the obtained iterative solution was extremely 
accurate, up to 10~ 15 , it was numerically confirmed the statement earlier proved and verified in Ref. || that 
once the quasilincar iteration sequence starts to converge, it will always continue to do so unlike the perturbation 
or 1/N expansion series, which are asymptotically divergent. The quasilinearization method therefore always 
yields the required precision once a successful initial guess generates convergence after a few steps. 

In view of all this, the quasilinearization method appears to be extremely useful in quantum mechanics and in 
many cases more advantageous than the perturbation theory or its different modifications, like expansion in inverse 
powers of the coupling constant, the 1/N expansion, etc. Though in this work only central potentials and one 
dimensional double well potential were considered and thus only differential equations in one variable were treated, 
the quailinearization method is able to solve the systems of nonlinear ordinary and partial differential equations in N 
variables and could therefore be applicable to the solution of the Schrodinger equation with the noncentral potentials 
or to the TV-body Schrodinger in 37V — 3 dimensions which will be subject of the future research. 
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TABLE I: QLM and exact binding energies E for different potentials. E are taken from citations in Ref. |l4j]. E p and N p are the 
energies by the 1/N perturbation method of Ref. [[HJ and the corresponding ranges of N where the 1/N expansion converges; 
a finite range means that the expansion diverges for larger N; stable digits are given only, n is the principal quantum number 
of the state. The uncertainty in last digit is in brackets where necessary for presentation, m denotes the (reduced) mass of the 
particle. 
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FIG. 1: Convergence of the logarithm of the absolute value of the difference of two successive QLM iterations Uk(r) for all r 
with the iteration index k for the ground state of the linear potential V = 2 7//2 r, m = 1. Here u(r) = arctan(— kxM/x'M) 
and k — \/2mE. The matching point is at r = 4. 




FIG. 2: As in Fig. |l|, but the convergence with respect to the exact solution obtained by solving the differential equation for 
0, Eq. (0. 
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FIG. 3: Convergence of the QLM iterations with the iteration index k for the solution of Fig. |l| Only a subset of iterations 
is presented such that the highest ones are not distinguishable from the exact solution. 




FIG. 5: As in Fig. |, but for the state of Fig. §. 




FIG. 8: As in Fig. |, but for the state of Fig. @. 



FIG. 10: As in Fig. [j], but for the ground (symmetric) state in the double-well potential V = (r 2 - R 2 ) 2 /(8R 2 ), R = 4 and 
m = 1/2. 




FIG. 11: As in Fig. Q, but the potential and state of Fig. hoi 
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FIG. 13: As in Fig. [l], but for the first excited (antisymmetric) state in the double-well potential V = (r 2 — R 2 ) 2 /(8R 2 ), R = 4 
and m = 1 /2. 




FIG. 14: As in Fig. Q, but the potential and state of Fig. [131 



